<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Example 6.4: Regressor selection problem</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/cvxbook/Ch06_approx_fitting/html/regressor_cvx.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Example 6.4: Regressor selection problem</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.3.1, Figure 6.7</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX Argyris Zymnis - 10/2005</span>
<span class="comment">%</span>
<span class="comment">% Solves</span>
<span class="comment">%        minimize   ||A*x-b||_2</span>
<span class="comment">%        subject to card(x) &lt;= k</span>
<span class="comment">%</span>
<span class="comment">% where card(x) denotes the number of nonzero elements in x,</span>
<span class="comment">% by first solving (for some value of alpha close to ||x_ln||_1)</span>
<span class="comment">%        minimize   ||A*x-b||_2</span>
<span class="comment">%        subject to ||x||_1 &lt;= alpha</span>
<span class="comment">%</span>
<span class="comment">% and iteratively decreasing alpha so as to get card(x) = k</span>
<span class="comment">% The sparsity pattern is then fixed in A and b and</span>
<span class="comment">%        minimize   ||A*x-b||_2</span>
<span class="comment">%</span>
<span class="comment">% is solved</span>

rand(<span class="string">'state'</span>,0);

m = 10;
n = 20;

A = randn(m,n);
b = A*[randn(round(m/2),1); zeros(n-round(m/2),1)];
b = b + 0.1*norm(b)*randn(m,1);

<span class="keyword">if</span> (1) <span class="comment">%%%%%%%%%%%%</span>

<span class="comment">% tradeoff curve for heuristic</span>
<span class="comment">%</span>
<span class="comment">% min.  ||Ax-b||_2</span>
<span class="comment">% s.t.  ||x||_1 &lt;= alpha</span>

residuals_heur = [norm(b)];
xln = A'*((A*A')\b);
lnorm = norm(xln,1);
nopts = 100;
alphas = linspace(0,lnorm,nopts);
residuals_heur = [norm(b)];
card_heur = [0];


<span class="keyword">for</span> k=2:(nopts-1)
  alpha = alphas(k);

  cvx_begin <span class="string">quiet</span>
    variable <span class="string">x(n)</span>
    minimize(norm(A*x-b))
    subject <span class="string">to</span>
        norm(x,1) &lt;= alpha;
  cvx_end

  x(find(abs(x) &lt; 1e-3*max(abs(x)))) = 0;
  ind = find(abs(x));
  sparsity = length(ind);
  fprintf(1,<span class="string">'Current sparsity pattern k = %d \n'</span>,sparsity);
  x = zeros(n,1);  x(ind) = A(:,ind)\b;
  card_heur = [card_heur, sparsity];
  residuals_heur = [residuals_heur, norm(A*x-b)];
<span class="keyword">end</span>;

obj1 = norm(b)
obj2 = [0];

i=1;
<span class="keyword">for</span> k=1:m-1
  <span class="keyword">if</span> ~isempty(find(card_heur == k))
     obj2(i+1) = k;
     obj1(i+1) = min(residuals_heur(find(card_heur ==k)));
     i=i+1;
  <span class="keyword">end</span>;
<span class="keyword">end</span>;
obj2(i) = m;  obj1(i) = 0;

<span class="keyword">end</span>; <span class="comment">%%%%%%%%%%%%%%%%%%%</span>


<span class="comment">% globally optimal tradeoff</span>


<span class="keyword">if</span> (1) <span class="comment">%%%%%%%%%%%%%</span>

bestx = zeros(n,m);
bestres = zeros(1,m);

<span class="keyword">for</span> k=1:m-1
  k
  <span class="comment">% enumerate sparsity patterns with exactly k nonzeros</span>
  bestres(k) = Inf;
  ind = 1:k
  nocases = 1;
  done = 0;
  <span class="keyword">while</span> ~done
     done = 1;
     <span class="keyword">for</span> i=0:k-1
       <span class="keyword">if</span> (ind(k-i) &lt; n-i),
          ind(k-i:k) = ind(k-i)+[1:i+1];
          done = 0;
          <span class="keyword">break</span>;
       <span class="keyword">end</span>;
     <span class="keyword">end</span>;
     <span class="keyword">if</span> done, <span class="keyword">break</span>; <span class="keyword">end</span>;
     x = zeros(n,1);
     x(ind) = A(:,ind)\b;
     <span class="keyword">if</span> (norm(A*x-b) &lt; bestres(k)),
        bestres(k) = norm(A*x-b);
        bestx(:,k) = x;
     <span class="keyword">end</span>;
     nocases = nocases + 1;
  <span class="keyword">end</span>;
  nocases
  factorial(n)/(factorial(n-k)*factorial(k))
<span class="keyword">end</span>;

x = A\b;
bestres(m) = norm(A*x-b);
bestres = [norm(b) bestres];

<span class="keyword">end</span>; <span class="comment">%%%%%%%%%</span>

figure
hold <span class="string">off</span>
obj1dbl =[];
obj2dbl =[];
<span class="keyword">for</span> i=1:length(obj1)-1
  obj1dbl = [obj1dbl, obj1(i), obj1(i)];
  obj2dbl = [obj2dbl, obj2(i), obj2(i+1)];
<span class="keyword">end</span>;
obj1dbl = [obj1dbl, obj1(length(obj1))];
obj2dbl = [obj2dbl, obj2(length(obj1))];

bestobj1 = bestres;
bestobj2 = [0:1:m];
bestobj1dbl =[];
bestobj2dbl =[];
<span class="keyword">for</span> i=1:length(bestobj1)-1
  bestobj1dbl = [bestobj1dbl, bestobj1(i), bestobj1(i)];
  bestobj2dbl = [bestobj2dbl, bestobj2(i), bestobj2(i+1)];
<span class="keyword">end</span>;
bestobj1dbl = [bestobj1dbl, bestobj1(length(bestobj1))];
bestobj2dbl = [bestobj2dbl, bestobj2(length(bestobj1))];

plot(obj1dbl,obj2dbl,<span class="string">'-'</span>, bestobj1dbl, bestobj2dbl,<span class="string">'--'</span>);
hold <span class="string">on</span>
plot(obj1,obj2,<span class="string">'o'</span>, bestobj1, bestobj2,<span class="string">'o'</span>);
axis([0 ceil(2*norm(b))/2 0 m+1])
xlabel(<span class="string">'x'</span>);
ylabel(<span class="string">'y'</span>);
hold <span class="string">off</span>

<span class="comment">%print -deps sparse_regressor_global_helv.eps</span>
<span class="comment">%save regressor_results</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 1 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 2 
Current sparsity pattern k = 3 
Current sparsity pattern k = 3 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 4 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 5 
Current sparsity pattern k = 6 
Current sparsity pattern k = 7 
Current sparsity pattern k = 8 
Current sparsity pattern k = 7 
Current sparsity pattern k = 7 
Current sparsity pattern k = 7 
Current sparsity pattern k = 9 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 10 
Current sparsity pattern k = 9 
Current sparsity pattern k = 9 
Current sparsity pattern k = 8 
Current sparsity pattern k = 8 
Current sparsity pattern k = 8 
Current sparsity pattern k = 8 
Current sparsity pattern k = 8 
Current sparsity pattern k = 9 
Current sparsity pattern k = 10 
Current sparsity pattern k = 11 
Current sparsity pattern k = 17 
Current sparsity pattern k = 17 
Current sparsity pattern k = 18 
Current sparsity pattern k = 18 
Current sparsity pattern k = 18 
Current sparsity pattern k = 18 
Current sparsity pattern k = 18 
Current sparsity pattern k = 18 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 19 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 
Current sparsity pattern k = 20 

obj1 =

    2.8954


k =

     1


ind =

     1


nocases =

    20


ans =

    20


k =

     2


ind =

     1     2


nocases =

   190


ans =

   190


k =

     3


ind =

     1     2     3


nocases =

        1140


ans =

        1140


k =

     4


ind =

     1     2     3     4


nocases =

        4845


ans =

        4845


k =

     5


ind =

     1     2     3     4     5


nocases =

       15504


ans =

       15504


k =

     6


ind =

     1     2     3     4     5     6


nocases =

       38760


ans =

       38760


k =

     7


ind =

     1     2     3     4     5     6     7


nocases =

       77520


ans =

       77520


k =

     8


ind =

     1     2     3     4     5     6     7     8


nocases =

      125970


ans =

      125970


k =

     9


ind =

     1     2     3     4     5     6     7     8     9


nocases =

      167960


ans =

      167960

</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="regressor_cvx__01.png" alt=""> 
</div>
</div>
</body>
</html>